﻿using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Drawing;
using System.Data;
using System.Text;
using System.Windows.Forms;
using SMcMaster;

namespace NumericalMethod
{
   public partial class LinearSystemForm : ContentPane
   {
      private StatusForm mStatus;
      private VParser FormulaParser = new VParser();//= new Parser();
      private List<List<NumericTextBox>> matrixBox = new List<List<NumericTextBox>>();
      private List<Label> equalLabel = new List<Label>();
      private List<NumericTextBox> equals = new List<NumericTextBox>();
      const int vx = 71, vy = 78, //offset matrix box
           sx = 41, sy = 20, //size textbox matrix
           gx = 2, gy = 3; //gap between matrix


      private int Ncount = 0;
      private int ActiveMethod = 0;

      private string[] method = new string[4] { 
           "Gaussian Elimination", 
           "Gaussian Elimination with Partial Pivoting", 
           "LU Decomposition",
           "Gauss-Seidel Method"
        };
      const string separator = "==================================";

      private string formatfPoint = "f4";

      public LinearSystemForm()
      {
         InitializeComponent();
         mStatus = new StatusForm();
      }

      private void Populate(int N)
      {
         if (N > 0)
         {
            for (int y = 0; y < Ncount + N; y++)
            {
               List<NumericTextBox> box;
               int x;
               if (matrixBox.Count > 0 && matrixBox.Count < Ncount || matrixBox.Count > y)
               {
                  box = matrixBox[y];
                  x = Ncount;
               }
               else
               {
                  box = new List<NumericTextBox>();
                  x = 0;
                  matrixBox.Add(box);
               }
               for (; x < Ncount + N; x++)
               {
                  NumericTextBox temp = new NumericTextBox();
                  temp.Size = new Size(sx, sy);
                  temp.Location = new Point(vx + (sx + gx) * x, vy + (sy + gy) * y);
                  Controls.Add(temp);
                  box.Add(temp);
               }

            }
         }
         else
         {
            for (int y = Ncount - 1; y >= 0; y--)
            {
               for (int z = Ncount - 1; z >= 0; z--)
               {
                  if (z + 1 == Ncount || y + 1 == Ncount)
                  {
                     Controls.Remove(matrixBox[y][z]);
                     matrixBox[y].Remove(matrixBox[y][z]);
                  }
               }
               if (matrixBox[y].Count == 0)
                  matrixBox.Remove(matrixBox[y]);
               //matrixBox.Remove(box);
            }
         }
         Ncount += N;
         addColumn();

         iterationLabel.Location = new Point(iterationLabel.Location.X, 20 + vx + (gy + sy) * Ncount);
         TabOrderManager tom = new TabOrderManager(this);
         tom.SetTabOrder(TabOrderManager.TabScheme.AcrossFirst);
      }

      public override void OnLoad()
      {
         Parent.Controls.Add(mStatus);
         methodComboBox.SelectedIndex = 0;
         mStatus.HidePanel(2);
         mStatus.Height = 400;
      }
      private void addColumn()
      {
         foreach (Label item in equalLabel)
            Controls.Remove(item);
         foreach (NumericTextBox item in equals)
            Controls.Remove(item);

         equalLabel.Clear();
         equals.Clear();
         for (int x = 0; x < Ncount; x++)
         {
            Label temp = new Label();
            temp.Size = new Size(20, sy);
            temp.Location = new Point(vx + (Ncount) * (sx + gx), vy + (x * (sy + gy)));
            temp.Text = "=";
            equalLabel.Add(temp);
            Controls.Add(temp);

            NumericTextBox tempz = new NumericTextBox();
            tempz.Size = new Size(sx, sy);
            tempz.Location = new Point(vx + 20 + (Ncount) * (sx + gx), vy + (x * (sy + gy)));
            equals.Add(tempz);
            Controls.Add(tempz);
         }
      }

      private void LinearSystemForm_Load(object sender, EventArgs e)
      {
         methodComboBox.Items.AddRange(method);
         Populate(2);
         methodComboBox.SelectedIndex = 0;
      }

      private void numericUpDown1_ValueChanged(object sender, EventArgs e)
      {
         Populate((int)numericUpDown1.Value - Ncount);
      }

      public float[,] MatrixValue
      {
         get
         {
            float[,] nnn = new float[Ncount, Ncount + 1];
            for (int x = 0; x < Ncount; x++)
            {
               for (int y = 0; y < Ncount; y++)
               {
                  nnn[y, x] = (float)matrixBox[y][x].DoubleValue;
               }
            }
            int index = 0;
            foreach (NumericTextBox txtBox in equals)
            {
               nnn[index++, Ncount] = (float)txtBox.DecimalValue;
            }
            return nnn;
         }
      }
      private void Calculate_Click(object sender, EventArgs e)
      {
         if (ActiveMethod == 0)
         {
            mStatus.AppendOutput("Gaussian Elimination");
            Gauss( Ncount, MatrixValue , false);
            mStatus.AppendOutput( separator);
         }
         else if ( ActiveMethod == 1 )
         {
            mStatus.AppendOutput("Gaussian Elimination with Partial Pivoting");
            Gauss( Ncount, MatrixValue, true );
            mStatus.AppendOutput( separator );
         }
         else if ( ActiveMethod == 2 )
         {
            mStatus.AppendOutput( "LU Decomposition" );
            int[] pivots = { 0 };
            //float[] ans;
            //matrixPrint( LU_Decomposition( MatrixValue, Ncount, Ncount, ref pivots ) );

            mStatus.AppendOutput( "Matrix form = " );
            matrixPrint( MatrixValue );

            DotNetMatrix.GeneralMatrix a = new DotNetMatrix.GeneralMatrix(
               MultiDimensionToJagged( MatrixValue, Ncount, Ncount ),
               Ncount,
               Ncount );
            mStatus.AppendOutput( "A = " );
            matrixPrint( GMToMultiDimension( a ) );

            DotNetMatrix.LUDecomposition LU = new DotNetMatrix.LUDecomposition( a );
            mStatus.AppendOutput( "L = " );
            matrixPrint( GMToMultiDimension( LU.L ) );
            mStatus.AppendOutput( "U = " );
            matrixPrint( GMToMultiDimension( LU.U ) );

            mStatus.AppendOutput( "Solution, UX = Y: " );
            DotNetMatrix.GeneralMatrix b = new DotNetMatrix.GeneralMatrix( getEqualColumn(), Ncount );
            matrixPrint( GMToMultiDimension( LU.Solve( b ) ) );
            //matrixPrint(ans);
            mStatus.AppendOutput( separator );
         }
         else if ( ActiveMethod == 3 )
         {
            float[] ans;
            mStatus.AppendOutput( "Gauss-Seidel Method" );
            try
            {
               SetFloatingPoint( (double) tolValue.DecimalValue );
               ans = Gauss_Seidel( MatrixValue, Ncount, (float) tolValue.DecimalValue, iterationValue.IntValue );
               answerPrint( ans );
            }
            catch ( Exception )
            {
               return;
            }
            mStatus.AppendOutput( separator );
         }
      }

      /// <summary>
      /// Set default floating point number
      /// for current calculation.
      /// </summary>
      /// <param name="err">Sample of double number</param>
      private void SetFloatingPoint(double err)
      {
         int fPoint = 4;//default is 4
         if (Math.Abs(err) < 1)
         {
            err = 1 / err;
            fPoint = 0;
            while (err > 1)
            {
               err /= 10;
               fPoint++;
            }
         }
         formatfPoint = "f" + fPoint;
      }
      static public double[][] MultiDimensionToJagged(float[,] a, int m, int n)
      {
         double[][] b = new double[m][];
         for (int i = 0; i < m; i++)
         {
            b[i] = new double[n];
            for (int j = 0; j < n; j++)
            {
               b[i][j] = a[i, j];
            }
         }
         return b;
      }
      public double[] getEqualColumn()
      {
         double[] e = new double[Ncount];
         for (int x = 0; x < Ncount; x++)
            e[x] = MatrixValue[x, Ncount];
         return e;
      }
      public float[,] GMToMultiDimension(DotNetMatrix.GeneralMatrix a)
      {
         double[][] ax = a.Array;
         float[,] b = new float[a.RowDimension, a.ColumnDimension];
         for (int i = 0; i < a.RowDimension; i++)
         {
            for (int j = 0; j < a.ColumnDimension; j++)
            {
               b[i, j] = (float)ax[i][j];
            }
         }
         return b;
      }
      private void matrixPrint(float[,] num)
      {
         string table;
         table = "\n<table border=\"1\" cellpadding=\"0\" cellspacing=\"0\">";
         for (int y = 0; y < num.GetLength(0); y++)
         {
            table += "<tr>";
            for (int x = 0; x < num.GetLength(1); x++)
            {
               table += ("<td>" + num[y, x] + "</td>");
            }
            table += "</tr>\n";

         }
         table += "</table>";
         mStatus.AppendHTML(table);
      }


      private void matrixPrint(float[] num)
      {
         string table;
         table = "\n<table border=\"1\" cellpadding=\"0\" cellspacing=\"0\">";
         for (int y = 0; y < num.GetLength(0); y++)
         {
            table += "<tr><td>" + num[y] + "</td></tr>\n";
         }
         table += "</table>";
         mStatus.AppendHTML(table);
      }

      private void methodComboBox_SelectedIndexChanged(object sender, EventArgs e)
      {
         ActiveMethod = methodComboBox.SelectedIndex;
         if (ActiveMethod == 3)
         {
            iterationLabel.Enabled = true;
            iterationValue.Enabled = true;
            tolLabel.Enabled = true;
            tolValue.Enabled = true;
            mStatus.HidePanel(0);
            chkMsgBoxInt msg = new chkMsgBoxInt();
            msg.ShowDialog();
         }
         else
         {
            iterationValue.Enabled = false;
            iterationLabel.Enabled = false;
            tolLabel.Enabled = false;
            tolValue.Enabled = false;
            mStatus.HidePanel(2);
         }
      }

      private void label2_LocationChanged(object sender, EventArgs e)
      {
         iterationValue.Location = new Point(iterationValue.Location.X, iterationLabel.Location.Y);
         tolLabel.Location = new Point(tolLabel.Location.X, iterationLabel.Location.Y);
         tolValue.Location = new Point(tolValue.Location.X, iterationLabel.Location.Y);
      }

      /// <summary>
      /// Solve system of N linear equations with N unknowns
      /// using Gaussian elimination with scaled partial pivoting
      /// First N rows and N+1 columns of A contain the system
      /// with right-hand sides of equations in column N+1
      /// err returns true if process fails; false if it is successful
      /// original contents of A are destroyed
      /// solution appears in column N
      /// </summary>
      /// <param name="N"></param>
      /// <param name="A"></param>
      /// <returns></returns>
      public float[] GaussianEllimation(int N, float[,] A)
      {
         int[] indx = new int[N * N];
         //N--;
         float[] scale = new float[N * N];
         float[] result = new float[N];
         float maxRatio;
         int maxIndx;
         int tmpIndx;
         float ratio;
         float sum;

         for (int i = 0; i < N; i++)
            indx[i] = i;	// index array initialization

         // determine scale factors

         for (int row = 0; row < N; row++)
         {
            scale[row] = Math.Abs(A[row, 0]);
            for (int col = 1; col < N; col++)
            {
               if (Math.Abs(A[row, col]) > scale[row])
                  scale[row] = Math.Abs(A[row, col]);
            }
         }

         // forward elimination

         for (int k = 0; k < N; k++)
         {
            // determine index of pivot row
            maxRatio = Math.Abs(A[indx[k], k]) / scale[indx[k]];
            maxIndx = k;
            for (int i = k + 1; i < N; i++)
            {
               if (Math.Abs(A[indx[i], k]) / scale[indx[i]] > maxRatio)
               {
                  maxRatio = Math.Abs(A[indx[i], k]) / scale[indx[i]];
                  maxIndx = i;
               }
            }
            if (maxRatio == 0) // no pivot available
            {
               throw new Exception("no pivot available");
            }
            tmpIndx = indx[k];
            indx[k] = indx[maxIndx];
            indx[maxIndx] = tmpIndx;

            // use pivot row to eliminate kth variable in "lower" rows
            for (int i = k + 1; i < N; i++)
            {
               ratio = -A[indx[i], k] / A[indx[k], k];
               for (int col = k; col <= N; col++)
               {
                  A[indx[i], col] += ratio * A[indx[k], col];
               }
            }
         }

         // back substitution

         for (int k = N - 1; k >= 0; k--)
         {
            sum = 0;
            for (int col = k + 1; col < N; col++)
            {
               sum += A[indx[k], col] * A[indx[col], N];
            }
            A[indx[k], N] = (A[indx[k], N] - sum) / A[indx[k], k];
         }


         for (int k = 0; k < N; k++)
            result[k] = A[indx[k], N];
         return result;
      }
      public void Gauss(int N, float[,] A, bool withPartialPivoting)
      {
          int[] indx = new int[N * N];
          float factor = 0;
          float[,] Disp = A;
          float[] result = new float[N];
          int nRow = N - 1;
          int tmpIndx = 0;
          int maxIndx = 0;

          matrixPrint(MatrixValue);
          
          for (int k = 0; k < N - 1; k++)
          {
              maxIndx = k;
              tmpIndx = indx[k];
              indx[k] = indx[maxIndx];
              indx[maxIndx] = tmpIndx;


              for (int i = k + 1; i < N + 1; i++)
              {
                  if (Disp[nRow, 0] == 0 && Disp[nRow, 1] == 0)
                      break;
                  if (i >= N)
                  {
                      if (Disp[nRow - 1, nRow - 1] == 0)
                      {
                          if (!withPartialPivoting)
                          {
                              mStatus.AppendOutput(" Please use Gauss Elimianation to avoid divided by zero");
                              return;
                          }

                      }

                      if (withPartialPivoting)
                      {
                          // Switched so that the largest element is the pivot element
                          //float[,] cDisp = Disp;
                          float[] lRow = new float[N + 1];
                          for (int p = 0; p < N + 1; p++)
                          {
                              lRow[p] = Disp[nRow - 1, p];
                          }
                          for (int p = 0; p < N + 1; p++)
                          {
                              Disp[nRow - 1, p] = Disp[nRow, p];
                          }
                          for (int p = 0; p < N + 1; p++)
                          {
                              Disp[nRow, p] = lRow[p];
                          }
                          mStatus.AppendOutput("Switch Row..");
                          matrixPrint(Disp);
                      }

                      factor = Disp[nRow, nRow - 1] / Disp[nRow - 1, nRow - 1];
                  }
                  else
                  {
                      if (A[k, k] == 0)
                      {
                          if (!withPartialPivoting)
                          {
                              mStatus.AppendOutput(" Please use Gauss Elimianation to avoid divided by zero");
                              return;
                          }
                          // Switched so that the largest element is the pivot element
                          float[,] cA = A;
                          for (int p = 0; p < N; p++)
                          {
                              A[i, p] = cA[k, p];
                          }
                      }


                      factor = A[i, k] / A[k, k];
                  }

                  if (factor > 0)
                  {
                      mStatus.AppendOutput("M :-" + factor.ToString());
                  }
                  else
                  {
                      mStatus.AppendOutput("M :" + factor.ToString());
                  }

                  //if (factor == 0)
                  //{
                  //    mStatus.AppendOutput(" Please use Gauss Elimianation to avoid divided by zero");
                  //    return;
                  //}

                  for (int j = k; j < N + 1; j++)
                  {

                      if (i >= N)
                      {
                          Disp[nRow, j] = Math.Abs(Disp[nRow, j] - factor * (Disp[nRow - 1, j])) < float.Parse("0.000001")
                              ? 0 : Disp[nRow, j] - factor * (Disp[nRow - 1, j]);
                      }
                      else
                      {
                          Disp[i, j] = Math.Abs(A[i, j] - factor * (A[k, j])) < float.Parse("0.000001") ? 0
                              : A[i, j] - factor * (A[k, j]);
                      }

                  }
                  matrixPrint(Disp);
              }
          }
          mStatus.AppendOutput("Result of Back Substitution:");
          Application.DoEvents();
          float sum = 0;
          for (int i = nRow; i >= 0; i--)
          {
              sum = Disp[i, nRow + 1];
              for (int j = i + 1; j < N; j++)
              {
                  sum = sum - (Disp[i, j] * result[j]);
              }
              result[i] = sum / Disp[i, i];
          }
          //for (int k = 0; k < N; k++)
          //    result[k] = Disp[indx[k], N];
          //matrixPrint(result);
          int r = 1;
          foreach (float i in result)
          {

              mStatus.AppendOutput("X" + r.ToString() + ":" + i.ToString());
              r++;
          }
      }
      public float[] Gauss_Seidel(float[,] Mtrx, int no, float tol, int iteration)
      {
         float[] X = new float[no]; // Matrix of values of X
         float[] oldX = new float[no];
         float error = float.MaxValue;
         List<string> columnList = new List<string>();
         
         //build column list
         columnList.Add("n");
         for (int i = 0; i < no; i++)
         {
            columnList.Add("x"+i);
         }
         columnList.Add("Error");
         mStatus.SetColumn(columnList.ToArray());

         for (int i = 0; i < no; i++)
         {
            Mtrx[i, no] = Mtrx[i, no] / Mtrx[i, i];
            for (int j = 0; j < no; j++)
            {
               if (i != j)
                  Mtrx[i, j] = Mtrx[i, j] / Mtrx[i, i];
            }
         }

         for (int k = 1; k <= iteration && error!=0 && error > tol; k++) // Applying Tolerance Limit
         {
            error = 0;
            Console.Write("ITE {0}\t", k);
            for (int i = 0; i < no; i++)
            {
               X[i] = Mtrx[i, no];

               for (int j = 0; j < no; j++)
               {
                  if (i != j)
                     X[i] = X[i] - Mtrx[i, j] * X[j];
               }
            }
            for (int i = 0; i < no; i++)
            {
               float newerror;
               Console.Write("X{0}={1:F3} \t", i, X[i]);
               newerror = Math.Abs(X[i]-oldX[i]);
               if (newerror > error)
                  error = newerror;
            }
            X.CopyTo(oldX,0);
            Console.WriteLine();

            //write row;
            columnList.Clear();//reuse back(bad programmer) :p
            columnList.Add(k.ToString());
            for (int i = 0; i < no; i++)
            {
               columnList.Add(X[i].ToString(formatfPoint));
            }
            columnList.Add(error.ToString(formatfPoint));
            mStatus.InsertRow(columnList.ToArray());
         }

         return X;
      }

      private void answerPrint(float[] num)
      {
         mStatus.AppendOutput("Answer: ");
         for (int y = 0; y < num.GetLength(0); y++)
         {
            mStatus.AppendOutput( "x" + (y + 1)+" = "+num[y]);
         }
         mStatus.AppendOutput(separator);
      }
      private void answerPrint( float[,] num )
      {
         mStatus.AppendOutput( "Answer: " );
         for ( int y = 0; y < num.GetLength( 1 ); y++ )
         {
            mStatus.AppendOutput( "x" + ( y + 1 ) + " = " + num[0, y ] );
         }
         mStatus.AppendOutput( separator );
      }

      /// <summary>
      /// LU decompostion.
      /// </summary>
      /// <param name="a"></param>
      /// <param name="m"></param>
      /// <param name="n"></param>
      /// <param name="pivots"></param>
      public static float[,] LU_Decomposition(float[,] a,
       int m,
       int n,
       ref int[] pivots)
      {
         int i = 0;
         int j = 0;
         int jp = 0;
         float[] t1 = new float[0];
         float s = 0;
         int i_ = 0;

         pivots = new int[Math.Min(m - 1, n - 1) + 1];
         t1 = new float[Math.Max(m - 1, n - 1) + 1];
         System.Diagnostics.Debug.Assert(m >= 0 & n >= 0,
            "Error in LUDecomposition: incorrect function arguments");

         //
         // Quick return if possible
         //
         if (m == 0 | n == 0)
         {
            return a;
         }
         for (j = 0; j <= Math.Min(m - 1, n - 1); j++)
         {

            //
            // Find pivot and test for singularity.
            //
            jp = j;
            for (i = j + 1; i <= m - 1; i++)
            {
               if (Math.Abs(a[i, j]) > Math.Abs(a[jp, j]))
               {
                  jp = i;
               }
            }
            pivots[j] = jp;
            if (a[jp, j] != 0)
            {

               //
               //Apply the interchange to rows
               //
               if (jp != j)
               {
                  for (i_ = 0; i_ <= n - 1; i_++)
                  {
                     t1[i_] = a[j, i_];
                  }
                  for (i_ = 0; i_ <= n - 1; i_++)
                  {
                     a[j, i_] = a[jp, i_];
                  }
                  for (i_ = 0; i_ <= n - 1; i_++)
                  {
                     a[jp, i_] = t1[i_];
                  }
               }

               //
               //Compute elements J+1:M of J-th column.
               //
               if (j < m)
               {
                  jp = j + 1;
                  s = 1 / a[j, j];
                  for (i_ = jp; i_ <= m - 1; i_++)
                  {
                     a[i_, j] = s * a[i_, j];
                  }
               }
            }
            if (j < Math.Min(m, n) - 1)
            {

               //
               //Update trailing submatrix.
               //
               jp = j + 1;
               for (i = j + 1; i <= m - 1; i++)
               {
                  s = a[i, j];
                  for (i_ = jp; i_ <= n - 1; i_++)
                  {
                     a[i, i_] = a[i, i_] - s * a[j, i_];
                  }
               }
            }
         }
         return a;
      }

   }
}
